function [elas_p, elas_p_mkt, major_unmerge] = elas_demand(root_dir, D, N, s1_predicted_y,...
    pr_y, N_y, resource_y1, resource_y2, ...
    alpha, variety_u, B1_c, B2_c, K, crra, rho, sigma_c1, sigma_c2, delta, beta_c, T1, T2)

% unmerging necessary and the relevant code is in supply estimation folder
cd(fullfile(root_dir, 'matlab/supply_estimation'));

% unmerge fringes
D_major = D.major;
D_nF = D.n_insurers(D_major==0); % N_st x 1, contains # of fringes in a given mkt

% unmerge market share variables
arg_share=1;
s1_unmerge = unmerge_fringe(D.s1, D_major, arg_share, D_nF, N.st); % NN x 1
s1_predicted_y = unmerge_fringe(s1_predicted_y', D_major, arg_share, D_nF, N.st)'; % N_y x NN

% unmerge other variables
arg_share=0;
p1 = unmerge_fringe(D.p1, D_major, arg_share, D_nF, N.st); % NN x 1
pr_s = unmerge_fringe(D.pr_s, D_major, arg_share, D_nF, N.st); % NN x K
major_unmerge = unmerge_fringe(D.major, D_major, arg_share, D_nF, N.st); % NN x 1
mkt = unmerge_fringe(D.mkt, D_major, arg_share, D_nF, N.st); % NN x 1

% store elasticity for each insurer in a given market
elas_p = zeros(length(p1),1);

for jj=1:length(p1)
    
    % set future rates to initial rates
    p2_input_jj = p1(jj)*ones(1, K);
    
    c1 = resource_y1-p1(jj); % N_y x 1
    c2 = repmat(resource_y2,1,K)-repmat(p2_input_jj,N_y,1); % N_y x K 
    
    % compute the (average) expected probability of the consumer staying w/ firm j
    agg_stay = pr_s(jj,:)*(1-delta)'; 
    agg_stay_sum=sum(agg_stay);
    
    D_v_p1 = (-B1_c*alpha*D_u_c(c1, rho, crra)...
        - ((beta_c^T1)*B2_c*alpha)*D_u_c(c2(:,1), rho, crra).*repmat(agg_stay_sum, N_y,1))/sigma_c1;
    
    D_s1_ind_p1 = repmat(p1(jj)', N_y, 1) .* D_v_p1.*(1-s1_predicted_y(:,jj)); % N_y x 1
    D_s1_p1 = sum(D_s1_ind_p1.*pr_y)'; % scalar
    
    elas_p(jj) =  D_s1_p1;
end

% compute weighted mean for each market using market shares
elas_p_mkt = zeros(N.st, 1);
for mm=1:N.st
     iind = find(mkt==mm); 
     elas_p_mm = elas_p(iind);
     s1_mm = s1_unmerge(iind);
     n_s1_mm = s1_mm/sum(s1_mm); % normalize market shares to be used as weights
     elas_p_mkt(mm) = sum(elas_p_mm.*n_s1_mm);
end
